2021: Gulf of Maine’s Warmest Year on Record

Characteristics of a record setting year for ocean temperatures

hide

# Set ggplot theme for figures
theme_set(theme_bw() + 
            theme(
              # Titles
              plot.title = element_text(hjust = 0, face = "bold", size = 14),
              plot.subtitle = element_text(size = 9),
              plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40"),
              legend.title = element_text(size = 9),
              legend.text = element_text(size = 7.5),
              # Axes
              axis.line.y = element_line(color = "black"),
              axis.ticks.y = element_line(), 
              axis.line.x = element_line(color = "black"),
              axis.ticks.x = element_line(), 
              axis.text = element_text(size = 11),
              rect = element_rect(fill = "transparent", color = "black"),
              # Facets
              strip.text = element_text(color = "white", 
                                        face = "bold",
                                        size = 11),
              strip.background = element_rect(
                color = "#00736D", 
                fill = "#00736D", 
                size = 1, 
                linetype="solid"))
          )



# Building a GMRI theme based on Wall street Journal and NYTimes theme
# base settings from {ggthemes}
theme_gmri <- function(base_size = 10, 
                       bg_color = "lightblue", 
                       base_family = "sans", 
                       title_family = "sans",
                       facet_color = "teal") {
    # Color from gmRi palette, sets background color
    #colorhex <- gmRi::gmri_cols()[bg_color]
    facet_hex <- gmri_cols()[facet_color]
    
    # Set up theme
    theme_foundation(
      base_size = base_size, 
      base_family = base_family) + 
        theme(
          
          # Major Elements
          line = element_line(linetype = 1, colour = "black"), 
          rect = element_rect(fill = "transparent", 
                              linetype = 0, 
                              colour = NA), 
          text = element_text(colour = "black"), 
          title = element_text(family = title_family, size = 12), 
          
          # Axis elements
          axis.text.x = element_text(colour = NULL), 
          axis.text.y = element_text(colour = NULL), 
          axis.ticks = element_line(colour = NULL), 
          axis.ticks.y = element_blank(), 
          axis.ticks.x = element_line(colour = NULL), 
          axis.line = element_line(), 
          axis.line.y = element_blank(), 
          axis.text = element_text(size = 11),
          
          # Legend Elements
          legend.background = element_rect(), 
          legend.position = "top", 
          legend.direction = "horizontal", 
          legend.box = "vertical", 
          legend.title = element_text(size = 9),
          legend.text = element_text(size = 7.5),
          
          # Panel/Grid Setup
          panel.grid = element_line(colour = NULL, linetype = 3, color = "gray70"), 
          panel.grid.major = element_line(colour = "black"), 
          panel.grid.major.x = element_blank(), 
          panel.grid.minor = element_blank(), 
          
          # Title and Caption Details
          plot.title = element_text(hjust = 0, face = "bold", size = 14),
          plot.subtitle = element_text(size = 9),
          plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
          plot.margin = unit(c(1, 1, 1, 1), "lines"), 
          
          # Facet Details
          strip.text = element_text(color = "white", face = "bold", size = 11),
          strip.background = element_rect(
            color = "white", 
            fill = facet_hex, 
            size = 1, 
            linetype="solid"))
}



# Set theme up for maps
map_theme <- list(
  theme(
    panel.border       = element_rect(color = "black", fill = NA),
    plot.background    = element_rect(color = "transparent", fill = "transparent"),
    line               = element_blank(),
    axis.title.x       = element_blank(), # turn off titles
    axis.title.y       = element_blank(),
    #legend.position    = "bottom", 
    legend.title.align = 0.5))

2021 Sea Surface Temperature

An Exceptional Year

2021 was an exceptionally warm year for the Gulf of Maine. It had the warmest fall on record for the area, and the second warmest summer. Much of the year the region experienced what would be considered “heatwave conditions”. After compiling sea surface temperature data for the entire uear we can now conclude that 2021 was the hottest year on Record for the Gulf of Maine.

hide

# Load the timeseries
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)


# Clean up the data - add labels
region_timeseries <- region_timeseries %>% 
  distinct(time, .keep_all = T) %>% 
  mutate(time = as.Date(time),
         yr = year(time),
         month_num = month(time),
         month = month(time, label = T, abbr = T),
         week = lubridate::week(time),
         doy = yday(time),
         season = metR::season(month_num, lang = "en"),
         #Set up correct year for winters, they carry across into next year
         season_eng = case_when(
           season == "SON" ~ "Fall",
           season == "DJF" ~ "Winter",
           season == "MAM" ~ "Spring",
           season == "JJA" ~ "Summer"),
         season_yr = ifelse( (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
         )

# Get heatwave statuses for each day:
# Uses area weighted sst by default
region_hw <- pull_heatwave_events(
  temperature_timeseries = region_timeseries, 
  threshold = 90, 
  clim_ref_period = c("1982-01-01", "2011-12-31"))

# Format dates and seasons and heatwave outcomes
region_hw <- region_hw %>% 
  mutate(time = as.Date(time),
         yr = year(time),
         month_num = month(time),
         month = month(time, label = T, abbr = T),
         week = lubridate::week(time),
         doy = yday(time),
         season = metR::season(month_num, lang = "en"),
         season_eng = case_when(
           season == "SON" ~ "Fall",
           season == "DJF" ~ "Winter",
           season == "MAM" ~ "Spring",
           season == "JJA" ~ "Summer"),
         #Set up correct year for winters, they carry across into next year
         season_yr = ifelse(
           (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
         )

# Set up heatwave status outcomes for fancy table features:
region_hw <- region_hw %>% 
  mutate(hw_outcomes = case_when(
    mhw_event == TRUE ~ 1,
    mcs_event == TRUE ~ 0,
    TRUE ~ 0.5))

Departure from Historic Conditions

Thanks to some long-running temperature monitoring programs we are able to track monthly sea surface temperatures as far back as 1850. The figure below shows how the region’s temperatures have fluctuated as measured by the Extended Reconstructed Sea Surface Temperature (ERSST) record.

Once we hit 1982 the availability of a higher resolution temperature resource becomes available in the form of the Optimum Interpolation Sea Surface Temperature (OISST) record. The ERSST data is of a coarser resolution (monthly measurements at a 0.5 x 0.5 degree grid) and does not capture the warmer inshore dynamics (a bias to show colder temperature), but is within half of a degree of the more modern measurements.

hide

# # Pathto ERSST on box
# ersst_path <- cs_path("okn", "ersst")
# 
# # Load temps and anoms (1982-2011 climatology)
# ersst_temp <- stack(str_c(ersst_path, "ERSSTv5_merged.nc"))
# ersst_anom <- stack(str_c(ersst_path, "ERSSTv5_anom.nc"))


# Updated 2022-02-4
ersst_path <- cs_path("res", "ERSSTv5")
# ersst_temp <- stack(str_c(ersst_path, "sst.mnmean_2022_02_02.nc")) #netcdf
ersst_tl <- read_csv(str_c(ersst_path, "ERSSTv5_anom_apershing_gulf_of_maine.csv")) 

# format
ersst_tl <- ersst_tl %>% 
  rename(Date = time) %>% 
  mutate(yr = year(Date),
         sst_clim = sst - sst_anom,
         area_wtd_clim = area_wtd_sst - area_wtd_sst_anom)


# Get Yearly Average
ersst_yr <- ersst_tl %>% 
  group_by(yr) %>% 
  summarise(sst = mean(sst, na.rm = T), 
            sst_anom = mean(sst_anom, na.rm = T),
            .groups = "drop") %>% 
  mutate(yr = as.numeric(yr)) %>% 
  filter(yr <= 2021) 


# Same for OISST
oisst_yr <- region_hw %>% 
  filter(between(yr, 1982, 2021)) %>% 
  group_by(yr) %>% 
  summarise(sst = mean(sst), .groups = "drop")

# Make Splines
ersst_smooth <- as.data.frame(spline(ersst_yr$yr, ersst_yr$sst))
oisst_smooth <- as.data.frame(spline(oisst_yr$yr, oisst_yr$sst))


# Plot them
ggplot() +
  geom_line(data = ersst_smooth, aes(x, y, color = "ERSSTv5"),
            group = 1, alpha = 0.4, linetype = 1) +
  geom_point(data = ersst_yr, aes(yr, sst, color = "ERSSTv5"), 
             size = 1) +
 geom_line(data = oisst_smooth, aes(x, y, color = "OISSTv2"),
            group = 1, alpha = 0.4, linetype = 1) +
  geom_point(data = oisst_yr, aes(yr, sst, color = "OISSTv2"), 
             size = 1) +
  scale_color_gmri() +
  theme_gmri() +
  labs(x = "", y = "Sea Surface Temperature", 
       color = "Temperature Data Record",
       title = "Long-Term Temperature Record for the Gulf of Maine",
       subtitle = "Current temperatures above 150-year highs") +
  scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C")) +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
  facet_zoom(xlim = c(1982, 2021)) +
  theme(legend.position = "bottom")

Daily Temperatures and Heatwave Status

hide

 # Last Calendar Year
this_yr <- region_hw %>% 
  filter(year(time) == 2021) 

#min/max of climate and min/max of 2021
clim_range <- round(range(this_yr$seas), 2)
clim_diff <- round(diff(clim_range), 2)

temp_range <- round(range(this_yr$sst), 2)
temp_diff <- round(diff(temp_range), 2)

anom_range <- round(range(this_yr$sst_anom), 2)
anom_diff <- round(diff(anom_range), 2)

The Gulf of Maine has a regular seasonal cycle hitting its lowest temperatures in March, and experiencing its peak highs in August. The difference between these two seasons is typically 22.914\(^oF\) or 12.73\(^oC\). In 2021 the temperature range was 5.09-20.49\(^oC\) or 15.4\(^oC\), with temperature anomaly values ranging from 0.84-4.09\(^oC\).

The highest temperatures anomalies spiked in the last days of June, capping off an already warm start to the year. A very-warm August and steady, above-average temperatures from early September through the end of the year helped solidifify 2021 as the hottest year on record for the Gulf of Maine.

hide

# Set colors by name
color_vals <- c(
  "Sea Surface Temperature" = "royalblue",
  "Heatwave Event"          = "darkred",
  "MHW Threshold"           = "coral3",
  "Daily Climatology"       = "gray30")


# Set the label with degree symbol
ylab <- "Sea Surface Temperature"



# Plot the last 365 days
hw_temp_p <- ggplot(this_yr, aes(x = time)) +
    geom_segment(aes(x = time, xend = time, y = seas, yend = sst), 
                 color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), 
                 color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    # geom_textpath(aes(y = mhw_thresh,  color = "MHW Threshold"), label = "Heatwave Threshold", hjust = 0.8, lty = 3) +
    #geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    geom_textpath(aes(y = seas,  color = "Daily Climatology"), label = "Climatology Mean", hjust = 0.5, lty = 2) +
    scale_color_manual(values = color_vals) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = expansion(mult = c(0,0))) +
    scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C")) +
    theme(legend.title = element_blank(),
          legend.position = "none") +
    labs(x = "", 
         y = ylab)



# Plot the last 365 days - anomaly scale
linetype_key <- c(
  "Sea Surface Temperature Anomaly" = 1,
  "Heatwave Event"                  = 1,
  "MHW Threshold"                   = 3,
  "Daily Climatology"               = 2)

# Same Plot as Anomalies
hw_anom_p <- this_yr %>% 
  mutate(
    sst = sst,
    seas = seas,
    sst_anom = sst_anom,
    mhw_thresh = mhw_thresh,
    anom_thresh = mhw_thresh - seas,
    anom_hwe = hwe - seas) %>% 
ggplot(aes(x = time)) +
    geom_segment(aes(x = time, xend = time, 
                     y = 0, yend = sst_anom), 
                 color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, 
                       y = anom_thresh, yend = anom_hwe), 
                 color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst_anom, color = "Sea Surface Temperature Anomaly")) +
    geom_line(aes(y = anom_hwe, color = "Heatwave Event")) +
    geom_line(aes(y = anom_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = 0, color = "Daily Climatology"), lty = 2, size = .5) +
    # geom_textpath(aes(y = 0,  color = "Daily Climatology"), 
    #               label = "Climatology Mean", hjust = 0.5, lty = 2, show.legend = FALSE) +
    scale_color_manual(values = color_vals) +
    scale_linetype_manual(values = linetype_key, guide = "none") +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = expansion(mult = c(0,0))) +
    guides(color = guide_legend(override.aes = list(linetype = linetype_key), nrow = 2)) +
    theme(legend.title = element_blank(),
          legend.position = "top") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                                           labels =  number_format(suffix = " \u00b0F")),
                       labels = number_format(suffix = " \u00b0C")) +
    labs(x = "", 
         y = "Temperature Anomaly", 
         caption = paste0("Climate reference period :  1982-2011"))


# Show figure
hw_temp_p / hw_anom_p

2021 Heatwave Events

hide

# Get the Event numbers
events_21 <-  this_yr %>% 
  drop_na(mhw_event_no) %>% 
  distinct(mhw_event_no) %>% 
  pull(mhw_event_no)

events_num <-  length(events_21)


# Pull their data, including the days before and after
events_dat <- filter(region_hw, mhw_event_no %in% as.character(events_21),
                     yr == 2021)
# Days in HW state
days_hw <- nrow(events_dat)

# How long on average
avg_duration <- round(days_hw / events_num, digits = 0)

During 2021 there were 5 main marine heatwave events, including one that carried over from the previous year and one that has continued into 2022. The average duration for these events was 72 days (excluding days in 2020 & 2022).

Record Daily Temperatures

To add context on the severity of these sustained heatwaves, many of the individual daily temperatures throughout 2021 were themselves the hottest days on record.

hide

# Grab the hottest day for each day of the year
hottest_days <- region_hw %>% 
  group_by(doy) %>% 
  slice_max(order_by = sst_anom, n = 1) %>% 
  ungroup() 

# Number in 2021
hottest_2021 <- hottest_days %>% 
  count(yr) %>% 
  filter(yr == "2021") %>% 
  pull(n)

# hottest_days %>% 
#   count(yr) %>% arrange(desc(n))

Of the 365 days in a year, 2021 set a new record temperature on 167 or 45.75% of the year.

Long-Term Heatwave Seasonality

When looking at the last ~40 years of marine heatwave events here in the Gulf of Maine. We see that the frequency, duration, and intensity of heatwave events has increased. This year was no exception with 358 days meeting the criteria for heatwave status.

hide

# Set up the date within a year
base_date <- as.Date("2000-01-01")
grid_data <- region_hw %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date))

# Prep the legend title
guide_lab <- "Sea Surface Temperature Anomaly \u00b0C"

# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) * c(-1, 1)


# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
  
  # background box fill for missing dates
  geom_rect(xmin = base_date, xmax = base_date + 365, 
            ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5, 
            fill = "gray75", color = "transparent") +
  
  # tile for sst colors
  geom_tile(aes(fill = sst_anom)) +
  
  # points for heatwave events
  geom_point(data = filter(grid_data, mhw_event == TRUE),
             aes(x = flat_date, y = year), size = .25)  +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
  scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
  labs(x = "", 
       y = "",
       "\nClimate reference period : 1982-2011") +
  
  #scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", 
                       limit = limit) +
  
  #5 inches is default rmarkdown height for barheight
  guides("fill" = guide_colorbar(title = guide_lab, 
                                 title.position = "right", 
                                 title.hjust = 0.5,
                                 barheight = unit(4.8, "inches"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +  
  theme(legend.title = element_text(angle = 90))


# Assemble pieces
heatwave_heatmap

Temperature Anomaly Horizons

One way to think about the severity of these changes is to think about temperature horizons. A temperature horizon captures how long temperatures remain above certain thresholds. Each threshold is designated its own temperature, and in this way we can see how long within a year temperatures remained: 1 to 2 or as much as 4\(^oC\) above normal.

hide

# make slight modification to grid_data
hori_data <- grid_data %>% 
  mutate(year = factor(year),
         year = fct_rev(year))


# cut out outliers and set cutpoints
cutpoints <- hori_data %>% 
  mutate(
    outlier = between(
      sst_anom,
      quantile(sst_anom, 0.25, na.rm=T)-
        1.5*IQR(sst_anom, na.rm=T),
      quantile(sst_anom, 0.75, na.rm=T)+
        1.5*IQR(sst_anom, na.rm=T))) %>% 
  filter(outlier)

# Set origin
ori <- 0

# Manually pick cut points
sca <- seq(-4, 4, length.out = 9)

sca_bins <- str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "\u00b0C")
sca_labels <- rev(sca_bins)


# Function to plot one or more years
anom_horizon_plot <- function(hori_data, origin = ori, scale_cutpoints = sca, labels = sca_labels){
  
   ggplot(hori_data) +
  geom_horizon(aes(flat_date, 
                   sst_anom,
                   fill = ..Cutpoints..), 
               origin = ori, 
               horizonscale = sca) +
  scale_fill_hcl(palette = 'RdBu', reverse = T, labels = sca_labels) +
  facet_grid(year~.) +
  theme(
    #strip.background = element_rect(fill = "#36454F", color = "#36454F"),
    panel.grid = element_blank(),
    panel.spacing.y=unit(0.1, "lines"),
    strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.border = element_blank(),
    legend.position = "left"
    ) +
  scale_x_date(expand=c(0,0), 
               date_breaks = "1 month", 
               date_labels = "%b") +
  labs(x = '', 
       fill = "Temperature Anomaly \u00b0C") 
  
  
}


# Plot Horizon Plot Using All Years
anom_horizon_plot(hori_data = filter(hori_data, yr <= 2021), 
                  origin = ori, 
                  scale_cutpoints = sca) +
  labs(title = "Temperature Anomaly Horizons",
       subtitle = "How deviation from the norm has become the new-normal",
       caption = "Climatatological Reference Period: 1982-2011")
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Comparing the Gulf’s Two Hottest Years

If we pull the horizons of our two hottest years on record it makes it easier to contrast where either one experienced acute high temperature events, and where there were sustained periods of above average temperatures.

Early into 2021 it was apparent that the year was on-par with the previous title-holder for warmest year on record.

Could have more text here

hide

# Reset the outliers
# cut out outliers and set cutpoints
cutpoints_2 <- hori_data %>% 
  filter(year %in% c("2012", "2021")) %>% 
  mutate(
    outlier = between(
      sst_anom,
      quantile(sst_anom, 0.25, na.rm=T)-
        1.5*IQR(sst_anom, na.rm=T),
      quantile(sst_anom, 0.75, na.rm=T)+
        1.5*IQR(sst_anom, na.rm=T))) %>% 
  filter(outlier)



# 2012
hori_2012 <- hori_data %>% 
  filter(year == "2012") %>% 
  anom_horizon_plot(scale_cutpoints = cutpoints_2) +
  facet_wrap(~year, strip.position = "top")
  
# 2021
hori_2021 <- hori_data %>% filter(year == "2021") %>% 
  anom_horizon_plot(scale_cutpoints = cutpoints_2) +
  facet_wrap(~year, strip.position = "top") +
  theme(plot.caption = element_text(colour = "gray25", size = 6)) +
  labs(caption = "Climatatological Reference Period: 1982-2011")

# Put them together as one figure
comparison_horizons <- hori_2012 / hori_2021 + plot_layout(guides = "collect")

# Do a bunch of formatting
comparison_horizons &
  theme(
    #Legend
    legend.position = "top",
    legend.key.height = unit(.5, "lines"),
    legend.key.width = unit(3.75, "lines")) &
   guides(fill = guide_legend(title = "Temperature Anomaly",
                              title.hjust = 0.5, 
                              title.position = "top",
                              label.position = "bottom", 
                              nrow = 1, 
                              byrow = T, reverse = T)) &
  plot_annotation(title = "Comparing the Gulf of Maine's Two Hottest Years", 
                  subtitle = "Duration of temperature anomaly horizons, colored by their strength")
  grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Losing the Balance Between Hot and Cold

Another way to visualize the climate transition that we are observing is by looking at the fraction of each year spent in different temperature ranges. Under a steady climate we would expect over the long-term to spend similar amounts of the year experiencing relatively warm & cold temperatures. These periods would balance themselves out and we would on-average have experienced something similar to the long-term climate.

What we have been seeing in the Gulf of Maine recently has lost that sense of balance. Larger fractions of the year are shared by above average temperatures & cold spell events are becoming vanishingly rare.

hide

# Set Factor Order for Bins
bin_order <- c(
  "Less than -1", 
  "-1 to -0.5",
  "Within 0.5",
  "0.5 to 1", 
  "Greater than 1")

# Set Colors
bin_colors <- RColorBrewer::brewer.pal(n = length(bin_order), 
                                       name = "RdBu")
bin_colors <- rev(bin_colors)

# Count Days in Bins
day_types <- region_hw %>% 
  mutate(
    `Less than -1`   = if_else(sst_anom <= -1, TRUE, FALSE),
    `-1 to -0.5`     = if_else(between(sst_anom, -1, -0.5), TRUE, FALSE),
    `Within 0.5`     = if_else(between(sst_anom, -0.5, 0.5), TRUE, FALSE),
    `0.5 to 1`       = if_else(between(sst_anom, 0.5, 1), TRUE, FALSE),
    `Greater than 1` = if_else(sst_anom >= 1, TRUE, FALSE)) %>% 
  group_by(yr) %>% 
  summarise(
    `Less than -1`   = sum(`Less than -1`, na.rm = T),
     `-1 to -0.5`    = sum(`-1 to -0.5`, na.rm = T),
    `Within 0.5`     = sum(`Within 0.5`, na.rm = T),
    `0.5 to 1`       = sum(`0.5 to 1`, na.rm = T),
    `Greater than 1` = sum(`Greater than 1`, na.rm = T),
    .groups = "drop") %>% 
  pivot_longer(cols = "Less than -1":"Greater than 1", names_to = "Anomaly Strength", values_to = "day_total", values_drop_na = FALSE) %>% 
  mutate(day_total = if_else(is.na(day_total), 0L, day_total),
         `Anomaly Strength` = factor(`Anomaly Strength`, levels = bin_order))


# Plot the Streamplot
hw_streamplot <- ggplot(day_types, aes(yr, y = day_total, fill = `Anomaly Strength`)) +
  geom_stream(type = "proportion") +
  #geom_segment(y = .5, yend = 0.5, linetype = 3, color = "gray50", x = 1981, xend = 2021, size = 0.25) +
  scale_x_continuous(breaks = seq(1980,2030, by = 10), 
                     minor_breaks = seq(1975, 2030, by = 5), 
                     expand = c(.03,0)) +
  scale_y_continuous(labels = scales::percent_format(), expand = c(0,0)) + 
  scale_fill_manual(values = setNames(bin_colors, bin_order)) + 
  annotate(x = 2015.5, y = 0.25, geom = "text", color = "white",
           label = "More than 50%\nof Days >1\u00b0C\nAbove Average") +
  labs(y = "", x = "", fill = "", 
       title = "A Warmer Gulf of Maine",
       subtitle = "Fractions of Each Year Experienced as Hot/Cold Anomalies",
       caption = "1981 not a complete years in the data*") +
  theme_minimal(base_size = 10) +
  theme(
    axis.line.x = element_line(),
    axis.ticks.x = element_line(),
    axis.line.y = element_line(),
    axis.ticks.y = element_line(),
    legend.position = "top",
    legend.key.height = unit(.5, "lines"),
    legend.key.width = unit(3.75, "lines"),
    panel.grid.minor = element_blank(),
    panel.grid = element_line(size = 0.2),
    plot.margin = unit(c(.5,1,.3,.5), "cm"),
    plot.title.position = "plot",
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 9),
    plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 7.5),
    legend.margin = margin(c(7,0,0,0)),
    legend.justification = "center") +
   guides(fill = guide_legend(title = "Anomaly Strength \u00b0C",
                             title.hjust = 0.5, 
                             title.position = "left",
                             label.position = "bottom", 
                             nrow = 1, 
                             byrow = T))


# Show figure
hw_streamplot
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Characteristics of the Hottest Years

Multivariate figure comparing characteristics of both years: 2012 & 2021

hide

# Use all years so they can be ranked independently

# Metrics?
year_mets <- hori_data %>% 
  filter(year != "2022") %>% 
  group_by(year) %>% 
  summarise(
    total_days                         = n(),
    `Average Temperature`              = round(mean(sst, na.rm = T), 2),
    `Temperature Above Average`        = round(mean(sst_anom, na.rm = T), 2),
    `Cumulative Degrees Above Average` = round(sum(sst_anom, na.rm = T), 0),
    hw_days                            = sum(mhw_event, na.rm = T),
    hw_events                          = n_distinct(mhw_event_no),
    peak_temp                          = max(sst),
    peak_anom                          = max(sst_anom)) %>% 
  mutate(
    hw_day_pct = round((hw_days/total_days) * 100, 2),
    `Average Heatwave Length` = round(hw_days/hw_events,  1)
  ) %>% 
  pivot_longer(names_to = "Var", values_to = "Metric", cols = `Average Temperature`:`Average Heatwave Length`)


# Snipe out the top5 for each
year_met_ranks <- year_mets %>% 
  filter(Var %in% c("Average Temperature", "Temperature Above Average", "Average Heatwave Length", "Cumulative Degrees Above Average")) %>% 
  group_by(Var) %>% 
  slice_max(n = 5, order_by = Metric) %>% 
  ungroup() %>%
  mutate(
    yr_focus = ifelse(year == "2021", TRUE, FALSE),
    year = reorder_within(year, Metric, Var)) 

# Plot the rankings
year_met_ranks %>% 
  ggplot(aes(year, Metric, color = yr_focus)) +
    geom_col(aes(fill = yr_focus), show.legend = F) +
    geom_label(aes(label = Metric),  show.legend = FALSE) +
    facet_wrap(~Var, scales = "free", ncol = 2) +
    coord_flip() +
    scale_x_reordered() +
    scale_y_continuous(expand = expansion(mult = c(0, .2))) +
    labs(y = "",
         x = NULL,
         title = "How 2021 Stacks Up:",
         subtitle = "Ranking Characteristics of Acute and Prolonged Warming Events") +
  theme(plot.title = element_text(hjust = 0)) + 
  scale_color_gmri() +
  scale_fill_gmri() +
  theme(panel.border = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major. = element_line(linetype = 3, color = "gray50"))
grid::grid.raster(lab_logo, x = 0.84, y = 0.94, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

---

Mapping Record Temperatures

2021 was an exceptionally warm year for many parts of the World, not just the Gulf of Maine. By taking the average temperature for the year for all years, we can rank how this year compares to previous years in the OISST data. For many places around the world this was one of the top 3 warmest years on record.

hide

# Load the rankings done in python, do some checking
ranks_path <- str_c(oisst_path, 'temp_rankings/yrly_ranks_1982to2021.nc')

#  Pull this year
ranks_21 <- stack(ranks_path)[["X2021"]]

# Reclassify Ranks
# Ranks are all integers so we just need to be sure that they are captured in the bounds we want
reclass_df <- c(
  0, 1.1, -1,  # Coldest Year
  1.1, 2.1, -2,  # Second Coldest
  2.1, 3.1, -3,  # Third Coldest
  3.1, 4.1, -4,  # Fourth Coldest
  4.1, 5.1, -5,  # Fifth Coldest
  3.1, 35, NA, # Things to Mask
  35.1, 36, 5, # Fifth Warmest
  36.1, 37, 4, # Fourth Warmest
  37.1, 38, 3, # Third Warmest
  38.1, 39, 2, # Second Warmest
  39.1, 40, 1  # Warmest
)

# make reclassification matrix
reclass_m <- matrix(reclass_df, ncol = 3, byrow = TRUE)

# reclassify - masking middle values handled with reclassification
ranks_reclass <- reclassify(ranks_21, reclass_m)

# Make Stars for ggplot
ranks_st <- st_as_stars(rotate(ranks_reclass))

# Mutate to get levels as factors for legend
ranks_st <- ranks_st %>% 
  mutate(
    temp_rank = case_when(
      X2021 == -1 ~ "Coldest",
      X2021 == -2 ~ "2nd Coldest",
      X2021 == -3 ~ "3rd Coldest",
      X2021 == -4 ~ "4th Coldest",
      X2021 == -5 ~ "5th Coldest",
      X2021 == 5 ~ "5th Warmest",
      X2021 == 4 ~ "4th Warmest",
      X2021 == 3 ~ "3rd Warmest",
      X2021 == 2 ~ "2nd Warmest",
      X2021 == 1 ~ "Warmest",
      TRUE ~ "Not a Record Year"),
    temp_rank = fct_drop(temp_rank),
    temp_rank = factor(temp_rank, levels = c(
      "Warmest", "2nd Warmest", "3rd Warmest", "4th Warmest", "5th Warmest",
      "Not a Record Year",
      "5th Coldest", "4th Coldest", "3rd Coldest", "2nd Coldest", "Coldest"))) %>% 
  select(-X2021) 


# Make Map
ranks_map <- ggplot() +
  geom_stars(data = ranks_st) +
  geom_sf(data = world_sf, size = .25) +
  scale_fill_brewer(palette = "RdBu", na.value = "transparent") +
  scale_x_continuous(breaks = seq(-180, 180, 90), 
                     expand = expansion(mult = c(0,0))) +
  scale_y_continuous(breaks = seq(-90, 90, 30), 
                     expand = expansion(mult = c(0,0))) +
  labs(title = "Record Setting Temperatures of 2021",
       fill = "",
       caption = "Temperature rankings cover period of 1982-2021") +
  map_theme
  
ranks_map
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

hide

# Zoom in on Northeast US
new_england <- ne_states("united states of america", returnclass = "sf")
canada      <- ne_states("canada", returnclass = "sf")



ggplot() +
  geom_stars(data = ranks_st) +
  geom_sf(data = new_england, size = .25) +
  geom_sf(data = canada, size = .25) +
  scale_fill_brewer(palette = "RdBu", na.value = "transparent")+ 
  coord_sf(xlim = c(-80, -55), ylim = c(34, 49)) +
  labs(x = "", y = "", fill = "2021 Temperature Rank",
       title = "Record Warm Temperatures",
       subtitle = "From North Carolina to Newfoundland",
       caption = "Temperature rankings cover period of 1982-2021")  +
  map_theme
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

hide

# Mapping in orthographic projection:

# Function for making a bounding box
make_bbox <- function(xlim = c(-75, -60), ylim = c(32, 46), crs =  4326){
 # Build bbox to crop with
  crop_bbox <- st_bbox(
    c(xmin = xlim[1],
      ymin = ylim[1],
      xmax = xlim[2],
      ymax = ylim[2]),
   crs = st_crs(crs))
  return(crop_bbox)
}


#zoom_bbox in different CRS
make_bbox() %>% 
  st_transform(ortho)



# # Issues: The warping step  is having issues, try cropping first
# crop_grid_bydim <- function(in_ras, xlim = c(-80, -50), ylim = c(20, 60), crs = crs){
#   
#   
#   
#   st_crop(in_ras, crop_bbox)
# }
# 
# 
# # Warp it to orthographic projection
# warp_to_ortho <- function(in_grid_st, ortho_lat = 25, ortho_lon = -50){
#   # orthographic projection
#   lat <- ortho_lon
#   lon <- ortho_lon
#   ortho <- paste0('+proj=ortho +lat_0=', lat, ' +lon_0=', lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')
#   
#   # Build grid in the crs you wish to warp to
#   projection_grid <- in_grid_st %>% 
#     st_transform(crs = ortho) %>% 
#     st_bbox() %>%
#     st_as_stars()
#   
#   # Warp to projection grid of chosen CRS
#   region_warp_ras <- in_grid_st %>% 
#     st_warp(projection_grid) 
#   
#   return(region_warp_ras)
#   
# }
# 
# 
# # Crop it down to a manageable area
# ranks_cropped <- crop_grid_bydim(in_ras = st_as_stars(rotate(ranks_21)), crs = 4326)
# 
# # use warping function:
# ranks_ortho <- warp_to_ortho(in_grid_st = ranks_cropped)


# source:
# https://gist.github.com/rafapereirabr/26965dd851debad32ad2e659024ba451


# Orthographic projection
ortho_lat <- 25
ortho_lon <- -60
ortho <- paste0('+proj=ortho +lat_0=', ortho_lat, ' +lon_0=', ortho_lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')


# Can we just project the original?
ranks_ortho <- projectRaster(rotate(ranks_21), crs = ortho)
ortho_st <- st_as_stars(ranks_ortho) %>% 
  mutate(
    temp_rank = case_when(
      X2021 == 1 ~ "Coldest",
      X2021 == 2 ~ "2nd Coldest",
      X2021 == 3 ~ "3rd Coldest",
      X2021 == 4 ~ "4th Coldest",
      X2021 == 5 ~ "5th Coldest",
      X2021 == 36 ~ "5th Warmest",
      X2021 == 37 ~ "4th Warmest",
      X2021 == 38 ~ "3rd Warmest",
      X2021 == 39 ~ "2nd Warmest",
      X2021 == 40 ~ "Warmest",
      TRUE ~ "Not a Record Year"),
    temp_rank = fct_drop(temp_rank),
    temp_rank = factor(temp_rank, levels = c(
      "Warmest", "2nd Warmest", "3rd Warmest", "4th Warmest", "5th Warmest",
      "Not a Record Year",
      "5th Coldest", "4th Coldest", "3rd Coldest", "2nd Coldest", "Coldest")))


# Try to:
# Fix polygons so they don't get cut in ortho projection
world <- rnaturalearth::ne_countries(scale = 'small', returnclass = 'sf')
world_ortho <- st_cast(world, 'MULTILINESTRING') %>%
  st_cast('LINESTRING', do_split=TRUE) %>%
  mutate(npts = mapview::npts(geometry, by_feature = TRUE)) %>%
  st_cast('POLYGON') %>% 
  st_transform(ortho) %>% 
  filter( is.na(st_dimension(.)) == FALSE )



# and map
ggplot() +
  geom_stars(data = select(ortho_st, -X2021)) +
  scale_fill_brewer(palette = "RdBu",  na.value = "gray50") +
  coord_sf(crs = ortho, expand = F) +
  map_theme

hide

# Preserving polygons in ortho view:
# https://gist.github.com/fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1

# checking geometries https://r-spatial.org/r/2017/03/19/invalid.html#empty-geometries

# Define the orthographic projection
# Choose lat_0 with -90 <= lat_0 <= 90 and lon_0 with -180 <= lon_0 <= 180
lat <- 45
lon <- -65
ortho <- paste0('+proj=ortho +lat_0=', lat, ' +lon_0=', lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')


# Load the world outlines
world <- ne_states(returnclass = "sf")



# Define the polygon that will help you finding the "blade"
# to split what lies within and without your projection
circle <- st_point(x = c(0,0)) %>% 
                   st_buffer(dist = 6371000) %>% 
                   st_sfc(crs = ortho)
  
# Project this polygon in lat-lon
circle_longlat <- circle %>% st_transform(crs = 4326)

# circle_longlat cannot be used as it is
# You must decompose it into a string with ordered longitudes
# Then complete the polygon definition to cover the hemisphere
if(lat != 0) {
    circle_longlat <- st_boundary(circle_longlat)
    
    circle_coords <- st_coordinates(circle_longlat)[, c(1,2)]
    circle_coords <- circle_coords[order(circle_coords[, 1]),]
    circle_coords <- circle_coords[!duplicated(circle_coords),]
    
    # Rebuild line
    circle_longlat <- st_linestring(circle_coords) %>% st_sfc(crs = 4326)
    
    if(lat > 0) {
      rectangle <- list(rbind(circle_coords,
                              c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
                              c(X = 180, Y = 90),
                              c(X = -180, Y = 90),
                              c(X = -180, circle_coords[1, 'Y']),
                              circle_coords[1, c('X','Y')])) %>% 
        st_polygon() %>% st_sfc(crs = 4326)
    } else {
      rectangle <- list(rbind(circle_coords,
                              c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
                              c(X = 180, Y = -90),
                              c(X = -180, Y = -90),
                              c(X = -180, circle_coords[1, 'Y']),
                              circle_coords[1, c('X','Y')])) %>% 
        st_polygon() %>% st_sfc(crs = 4326)
    }
    
    circle_longlat <- st_union(st_make_valid(circle_longlat), st_make_valid(rectangle))
  }
  
# This visualization shows the visible emisphere in red
ggplot() +
    geom_sf(data = world) +
    geom_sf(data = circle_longlat, color = 'red', fill = 'red', alpha = 0.3)
  


# A small negative buffer is necessary to avoid polygons still disappearing in a few pathological cases
# I should not change the shapes too much
  visible <- st_intersection(st_make_valid(world),   
                             st_buffer(circle_longlat, -0.09)) %>%
    st_transform(crs = ortho)

# DISCLAIMER: This section is the outcome of trial-and-error and I don't claim it is the best approach 
# Resulting polygons are often broken and they need to be fixed
# Get reason why they're broken
broken_reason <- st_is_valid(visible, reason = TRUE)
  
# First fix NA's by decomposing them
# Remove them from visible for now
na_visible <- visible[is.na(broken_reason),]
visible <- visible[!is.na(broken_reason),]
  
# Open and close polygons
na_visible <- st_cast(na_visible, 'MULTILINESTRING') %>% 
    st_cast('LINESTRING', do_split=TRUE)
na_visible <- na_visible %>% 
  mutate(npts = mapview::npts(geometry, by_feature = TRUE))

# Exclude polygons with less than 4 points
na_visible <- na_visible %>%
    filter(npts >=4) %>%
    select(-npts) %>%
    st_cast('POLYGON')

# Fix other broken polygons
broken <- which(!st_is_valid(visible))
for(land in broken) {
    result = tryCatch({
    # visible[land,] <- st_buffer(visible[land,], 0) # Sometimes useful sometimes not
    visible[land,] <- st_make_valid(visible[land,]) %>% 
        st_collection_extract()  
    }, error = function(e) {
      visible[land,] <<- st_buffer(visible[land,], 0)
    })
}

# Bind together the two tables
visible <- rbind(visible, na_visible)

# Final plot
ggplot() +
  geom_sf(data = circle,
        #fill = 'aliceblue') + # if you like the color
        fill = NA) +
    geom_sf(data = st_collection_extract(visible)) +
    coord_sf(crs = ortho)

Comparing Regions

We are often asked when speaking to the degree of the Gulf of Maine is warming: > “Where is it warming faster?”

This question is challenging because the Gulf of Maine, and any area in the ocean has a unique size/shape/dynamic that makes them difficult to directly compare directly. To put the Gulf of Maine’s warming into perspective against similar ecologically and oceanographically relevant units we can compare it against the large marine ecosystems (LME) of the world.

hide

# Load the LME timeseries to calculate rates of warming 
lme_names <- get_region_names("lme") # names
lme_paths <- get_timeseries_paths("lme", mac_os = "mojave") # paths

# Data
lme_oisst <- map(lme_names, ~ read_csv(lme_paths[[.x]][["timeseries_path"]], col_types = cols(), guess_max = 1e6)) 
lme_oisst <- setNames(lme_oisst, lme_names)

# Add Gulf of Maine Timeseries to the List with LME's
lme_oisst[["gulf_of_maine"]] <- region_timeseries

# Make sure dates are formatted
lme_oisst <- map(lme_oisst, ~ mutate(.x, time = as.Date(time)))


####  Warming Rates  ####
lme_8221 <- map(lme_oisst, ~ mutate(.x, year = year(time)) %>% filter(year %in% c(1982:2021)))

# Get warming rates
lme_rates <- map_dfr(lme_8221, function(lme_data){
  # Uses area weighted SST
  lme_yearly <- group_by(lme_data, year) %>% 
    summarise(sst = mean(area_wtd_sst, na.rm = T))
  yrly_warming <- lm(sst ~ year, data = lme_yearly)
  mod_details <- broom::tidy(yrly_warming)
  yrly_rate <- mod_details[2,"estimate"]
  names(yrly_rate) <- "annual_warming_rate_C"
  yrly_rate
}, .id = "Region") %>% 
  arrange(desc(annual_warming_rate_C))


# Pull out the top 13
rates_ordered <- lme_rates %>% 
  mutate(`Warming Rank` = row_number(),
         `Warming Rate F` = as_fahrenheit(annual_warming_rate_C, data_type = "anomalies"),
         Region = str_replace_all(Region, "_", " "),
         Region = str_to_title(Region)) %>% 
  select(Region, `Warming Rank`, `Warming Rate C` = annual_warming_rate_C, `Warming Rate F`) 

# symbol for degree: \u00b0

hide

# Make table of top 10*
rates_ordered %>% 
  slice(1:10) %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  gt(rowname_col = "Region") %>% 
   tab_header(
    title = md("**Warming Rates of Large Marine Ecosystems**")) %>%
  tab_stubhead(label = "Region") %>% 
  tab_style(
    style = list(
      cell_fill(color = "gray90"),
      cell_text(style = "oblique")),
    locations = cells_body(
      rows = Region == "Gulf Of Maine")) %>% 
  tab_spanner(label = md("*Annual Temperature Change*"),
              columns = c(`Warming Rate C`, `Warming Rate F`)) %>% 
 
  cols_label(
      `Warming Rate F` = "\u00b0F",
      `Warming Rate C` = "\u00b0C") %>% 
  tab_source_note(source_note = md("**Data Source:** *NOAA OISSTv2 Daily Sea Surface Temperature Data.*")) %>%                    
  tab_source_note(source_note = md("**Notes:** *Temperatures adjusted for latitudinal distotion.*") ) %>%  
  tab_footnote(
    footnote = "Gulf of Maine not a true large marine ecosystem",
    locations = cells_stub(rows = which(rates_ordered$Region == "Gulf Of Maine")) 
  )

Warming Rates of Large Marine Ecosystems
Region	Warming Rank	Annual Temperature Change
		°C	°F
Baltic Sea	1	0.047	0.085
Gulf Of Maine1	2	0.044	0.080
Black Sea	3	0.044	0.079
Scotian Shelf	4	0.041	0.073
Iceland Shelf And Sea	5	0.038	0.069
Northeast Us Continental Shelf	6	0.038	0.068
North Sea	7	0.035	0.063
Norwegian Sea	8	0.032	0.057
Sea Of Japan	9	0.031	0.056
Mediterranean Sea	10	0.030	0.055
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Notes: Temperatures adjusted for latitudinal distotion.
1 Gulf of Maine not a true large marine ecosystem

hide

# Netcdf of warming rates
rates_path <- str_c(oisst_path, "warming_rates/annual_warming_rates1982to2021.nc")
rates_nc <- stack(rates_path, varname = "annual_warming_rate")
ranks_nc <- stack(rates_path, varname = "rate_percentile")



# Get paths to the LME's and their timeseries
lme_shapes <- map(lme_names, ~  read_sf(lme_paths[[.x]][["shape_path"]]))  %>% 
  setNames(lme_names)

# Include Gulf of Maine as a shape
gom_shape <- read_sf(poly_path)
lme_shapes[["gulf_of_maine"]] <- gom_shape



# Masking Function to clip to the study area
mask_nc <- function(ras_obj, mask_shape){
  
  # Get stats from ranks
  m1 <- mask(rotate(ras_obj), mask_shape)
  m1 <- crop(m1, mask_shape)
  return(m1)
}



# Mask with each to get warming %
lme_rates_masked <- map(lme_shapes, ~ mask_nc(ras_obj = rates_nc, mask_shape = .x))
lme_perc_masked <- map(lme_shapes, ~ mask_nc(ras_obj = ranks_nc, mask_shape = .x))


# Put the percentiles and average rates in one table
if( all(names(lme_rates_masked) == names(lme_perc_masked)) ){
  percentile_table <- map2_dfr(lme_rates_masked, lme_perc_masked, function(rate_lme, perc_lme){
    rate_val <- cellStats(rate_lme, mean, na.rm = T)
    perc_val <- cellStats(perc_lme, mean, na.rm = T)
    df_out <- data.frame(
      "annual_rate" = rate_val,
      "warming_percentile" = perc_val)
    
    return(df_out)
    
  }, .id = "ecosystem")
}


# # Check it out:
# percentile_table %>% arrange(desc(warming_percentile)) %>% head()
# percentile_table %>% arrange(desc(annual_rate)) %>% head()


# Bizarre...

# Put LME's in one collection
lme_shape_tab <- bind_rows(lme_shapes)

# Join with the  rates/ranks
perc_tab <- percentile_table %>% 
  mutate(name = str_replace_all(ecosystem, "_", " "),
         name = str_to_title(name),
         name = str_replace(name, "Of", "of"))

# Add the percentile information, and map them
lme_sf_ranks <- left_join(lme_shape_tab, perc_tab, by = c("name")) 

ggplot() +
  geom_sf(data = lme_sf_ranks, aes(fill = warming_percentile), color = "transparent") +
  geom_sf(data = world_sf, size = 0.25) +
  scale_fill_distiller(palette = "RdBu", direction = -1, na.value = "transparent") +
  theme(legend.position = "bottom") +
  coord_sf(expand = FALSE) +
  map_theme +
  guides(fill = guide_colorbar(
    title.position = "top",
    title.hjust = 0.5, 
    barwidth = unit(4, "in"), 
    frame.colour = "black")) +
  labs(fill = "Warming Rate Percentile",
       title = "")

Biological Events

 

---

A work by Adam A. Kemberling

Akemberling@gmri.org